1 Introduction

Here, we will annotate the cells of the patient with id “365”.

2 Pre-processing

2.1 Load packages

library(Seurat)
library(SeuratWrappers)
library(harmony)
library(tidyverse)

2.2 Define parameters

# Paths
path_to_obj <- here::here("results/R_objects/5.seurat_clustered_365.rds")
path_to_obj_sampl_artifacts <- here::here("results/R_objects/cll_seurat_annotated.rds")
path_to_save <- here::here("results/R_objects/6.seurat_annotated_365.rds")
path_to_save_markers <- here::here("3-clustering_and_annotation/markers_clusters_365.rds")


# Colors
color_palette <- c("black", "gray", "red", "yellow", "violet", "green4",
                   "blue", "mediumorchid2", "coral2", "blueviolet",
                   "indianred4", "deepskyblue1", "dimgray", "deeppink1",
                   "green", "lightgray", "hotpink1")


# Source functions
source(here::here("bin/utils.R"))


# Thresholds
min_log2FC <- 0.3
alpha <- 0.001

2.3 Load data

seurat <- readRDS(path_to_obj)
seurat
## An object of class Seurat 
## 23326 features across 4685 samples within 1 assay 
## Active assay: RNA (23326 features, 2000 variable features)
##  3 dimensional reductions calculated: pca, harmony, umap
DimPlot(seurat, cols = color_palette)

seurat_bias <- readRDS(path_to_obj_sampl_artifacts)

3 Cell Cycle Score

As upregulation of cell cycle genes is a hallmark of Richter transformation, we will infer the cell cycle score and phase for each cell:

DimPlot(seurat, group.by = "time_point")

FeaturePlot(seurat, "CD3D")

seurat <- CellCycleScoring(
  seurat,
  s.features = cc.genes.updated.2019$s.genes,
  g2m.features = cc.genes.updated.2019$g2m.genes,
  set.ident = FALSE
)
DimPlot(seurat, group.by = "Phase")

umap_s_score <- FeaturePlot(seurat, features = "S.Score") +
  scale_color_viridis_c(option = "magma") +
  labs(title = "S Score") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12, face = "plain"),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.line = element_blank()
  )
umap_g2m_score <- FeaturePlot(seurat, features = "G2M.Score") +
  scale_color_viridis_c(option = "magma") +
  labs(title = "G2M Score") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12, face = "plain"),
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.line = element_blank()
  )
umap_cc_combined <- ggpubr::ggarrange(
  plotlist = list(umap_s_score, umap_g2m_score),
  nrow = 2,
  ncol = 1,
  common.legend = FALSE
)
umap_cc_combined

4 Find Markers

markers <- FindAllMarkers(seurat, only.pos = TRUE, logfc.threshold = min_log2FC)
markers <- markers %>%
  mutate(cluster = as.character(cluster)) %>%
  filter(p_val_adj < alpha) %>%
  arrange(cluster) %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE)
DT::datatable(markers, options = list(scrollX = TRUE))

5 Annotation

Important literature to annotate the cells:

Cluster Markers Annotation
0 CXCR4 CLL
1_0 CD27 CXCR4loCD27hi RT
1_1 MIR155HG MIR155HGhi RT
1_2 ENO1, TCL1A RT quiescent
1_3 PCNA, TOP2A RT proliferative

Annotate:

seurat$annotation_final <- factor(
  seurat$final_clusters,
  levels = c("0", "1_0", "1_1", "1_2", "1_3")
)
new_levels_365 <- c("CLL", "CXCR4loCD27hi RT", "MIR155HGhi RT",
                    "RT quiescent", "RT proliferative")
levels(seurat$annotation_final) <- new_levels_365
seurat$annotation_final <- factor(
  seurat$annotation_final,
  color_annotations$`365`$annotation_final
)
Idents(seurat) <- "annotation_final"


# Plot UMAP
(umap_annotation <- plot_annotation(
  seurat_obj = seurat,
  pt_size = 0.5,
  colors_reference = color_annotations,
  patient_id = "365",
  nothing = TRUE
))

6 Visualize markers

# UMAPs
genes_interest <- c("CXCR4", "CD27", "MIR155HG", "ENO1", "TCL1A",
                    "PCNA", "TOP2A")
(feature_plots <- purrr::map(genes_interest, function(x) {
  p <- FeaturePlot(seurat, x, pt.size = 0.5) +
    scale_color_viridis_c(option = "magma")
  p
}))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

# Dot plots
(dot_plot <- plot_dot_plot(
  seurat,
  goi = rev(genes_interest),
  colors_reference = color_annotations,
  patient_id = "365"
))

# Violin plots
(vln_plot_s <- plot_violin_plot(
  seurat,
  continuous_var = "S.Score",
  ylab = "S Phase Score",
  colors_reference = color_annotations,
  patient_id = "365"
))

(vln_plot_g2m <- plot_violin_plot(
  seurat,
  continuous_var = "G2M.Score",
  ylab = "G2M Phase Score",
  colors_reference = color_annotations,
  patient_id = "365"
))

7 Validation 1: CLL-specific clustering

Although we could recover the major clusters that we have observed in previous patients, it is intriguing that the clusters of MIR155HGhi and CXCR4loCD27hi were found in the RT time-point and not in the CLL. Thus, let us perform a CLL-specific clustering:

seurat_cll <- subset(seurat, time_point == "T2")
seurat_cll <- process_seurat(seurat_cll, dims = 1:20)
feature_plots2 <- purrr::map(c("CXCR4", "CD27", "MIR155HG", "S100A4"), function(x) {
  p <- FeaturePlot(seurat_cll, x, pt.size = 0.5) +
    scale_color_viridis_c(option = "magma")
  p
})
feature_plots2
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Indeed, we see that CXCR4, CD27 and MIR155HG are mutually exclusive, but we cannot distinguish 3 clear clusters

seurat_cll <- FindNeighbors(seurat_cll, dims = 1:20, reduction = "pca")
seurat_cll <- FindClusters(seurat_cll, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4119
## Number of edges: 134238
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6852
## Number of communities: 4
## Elapsed time: 0 seconds
DimPlot(seurat_cll)

8 Validation 2: sampling artifacts

Could this pattern be explained by a differential processing time between both samples (RT and CLL)?

seurat_bias <- subset(seurat_bias, temperature == "RT")
seurat_bias <- subset(seurat_bias, time == "24h")
DimPlot(seurat_bias)

seurat_bias <- SplitObject(seurat_bias, "cell_type")
seurat_bias <- seurat_bias[c("CLL 1472", "CLL 1892")]
seurat_bias <- purrr::map(seurat_bias, process_seurat, dims = 1:20)
DimPlot(seurat_bias$`CLL 1472`, group.by = "time")

FeaturePlot(seurat_bias$`CLL 1472`, c("CXCR4", "CD27", "MIR155HG"))

DimPlot(seurat_bias$`CLL 1892`, group.by = "time")

FeaturePlot(seurat_bias$`CLL 1892`, c("CXCR4", "CD27", "MIR155HG"))

9 Save

# Save Seurat object
saveRDS(seurat, path_to_save)


# Save markers
markers$annotation <- factor(markers$cluster)
levels(markers$annotation) <- new_levels_365
markers_list <- purrr::map(levels(markers$annotation), function(x) {
  df <- markers[markers$annotation == x, ]
  df <- df[, c(7, 1, 5, 2:4, 6, 8)]
  df
})
names(markers_list) <- levels(markers$annotation)
markers_list <- markers_list[color_annotations$`365`$annotation_final]
markers_final <- bind_rows(markers_list)
saveRDS(markers_list, path_to_save_markers)
saveRDS(
  markers_final,
  here::here("results/tables/markers/markers_annotated_clusters_patient_365.rds")
)
saveRDS(markers_list, path_to_save_markers)
openxlsx::write.xlsx(
  x = markers_list,
  file = "results/tables/markers/markers_annotated_clusters_patient_365.xlsx"
)

10 Session Info

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.6          purrr_0.3.4          readr_1.4.0          tidyr_1.1.3          tibble_3.1.2         ggplot2_3.3.3        tidyverse_1.3.1      harmony_1.0          Rcpp_1.0.6           SeuratWrappers_0.3.0 SeuratObject_4.0.2   Seurat_4.0.3         BiocStyle_2.18.1    
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1          backports_1.2.1       plyr_1.8.6            igraph_1.2.6          lazyeval_0.2.2        splines_4.0.4         crosstalk_1.1.1       listenv_0.8.0         scattermore_0.7       digest_0.6.27         htmltools_0.5.1.1     fansi_0.5.0           magrittr_2.0.1        tensor_1.5            cluster_2.1.1         ROCR_1.0-11           limma_3.46.0          openxlsx_4.2.3        remotes_2.4.0         globals_0.14.0        modelr_0.1.8          matrixStats_0.59.0    spatstat.sparse_2.0-0 colorspace_2.0-1      rvest_1.0.0           ggrepel_0.9.1         haven_2.4.1           xfun_0.23             crayon_1.4.1          jsonlite_1.7.2        spatstat.data_2.1-0   survival_3.2-10       zoo_1.8-9             glue_1.4.2            polyclip_1.10-0       gtable_0.3.0          leiden_0.3.8          car_3.0-10            future.apply_1.7.0    abind_1.4-5           scales_1.1.1          DBI_1.1.1             rstatix_0.7.0         miniUI_0.1.1.1        viridisLite_0.4.0     xtable_1.8-4          reticulate_1.20       spatstat.core_2.1-2   foreign_0.8-81        rsvd_1.0.5            DT_0.18               htmlwidgets_1.5.3     httr_1.4.2            RColorBrewer_1.1-2   
##  [55] ellipsis_0.3.2        ica_1.0-2             farver_2.1.0          pkgconfig_2.0.3       sass_0.4.0            uwot_0.1.10           dbplyr_2.1.1          deldir_0.2-10         here_1.0.1            utf8_1.2.1            labeling_0.4.2        tidyselect_1.1.1      rlang_0.4.11          reshape2_1.4.4        later_1.2.0           munsell_0.5.0         cellranger_1.1.0      tools_4.0.4           cli_2.5.0             generics_0.1.0        broom_0.7.7           ggridges_0.5.3        evaluate_0.14         fastmap_1.1.0         yaml_2.2.1            goftest_1.2-2         knitr_1.33            fs_1.5.0              fitdistrplus_1.1-5    zip_2.2.0             RANN_2.6.1            pbapply_1.4-3         future_1.21.0         nlme_3.1-152          mime_0.10             xml2_1.3.2            compiler_4.0.4        rstudioapi_0.13       curl_4.3.1            plotly_4.9.4          png_0.1-7             ggsignif_0.6.2        spatstat.utils_2.2-0  reprex_2.0.0          bslib_0.2.5.1         stringi_1.6.2         highr_0.9             RSpectra_0.16-0       lattice_0.20-41       Matrix_1.3-4          vctrs_0.3.8           pillar_1.6.1          lifecycle_1.0.0       BiocManager_1.30.15  
## [109] spatstat.geom_2.1-0   lmtest_0.9-38         jquerylib_0.1.4       RcppAnnoy_0.0.18      data.table_1.14.0     cowplot_1.1.1         irlba_2.3.3           httpuv_1.6.1          patchwork_1.1.1       R6_2.5.0              bookdown_0.22         promises_1.2.0.1      rio_0.5.26            KernSmooth_2.23-18    gridExtra_2.3         parallelly_1.26.0     codetools_0.2-18      MASS_7.3-53.1         assertthat_0.2.1      rprojroot_2.0.2       withr_2.4.2           sctransform_0.3.2     mgcv_1.8-36           parallel_4.0.4        hms_1.1.0             grid_4.0.4            rpart_4.1-15          rmarkdown_2.8         carData_3.0-4         Rtsne_0.15            ggpubr_0.4.0          shiny_1.6.0           lubridate_1.7.10